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Abstract 

This work deals with the integration of nonsmooth flexible multibody systems with impacts and dry friction. 
We develop a framework which improves a non-impulsive trajectory of state variables by impulsive correction 
after each time-step if necessary. This correction is automatic and is evaluated on the same kinematic level 
as the piecewise non-impulsive trajectory. The resulting overall mixed timestepping scheme is consistent with 
respect to impacts and friction as well as benefits from advantages of the base integration schemes used to 
calculate the approximation inside the time-step. Therefore, we compare the generalized-a method, the Bathe 
method and the ED-a method. 
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1 Introduction 

Nonsmooth mechanical systems are characterized by jumps in the system’s velocities or accelerations. Typical 
examples with impacts or dry friction can e.g. be found in automotive, railway or robotics applications l[mi2^ 
I2l[29l|34ll26l|47l. Hence, we deal with impulsive and non-impulsive periods. There are two types of integration 
methods to handle nonsmooth motion: event-driven schemes and timestepping schemes. On the one hand, 
event-driven schemes, i.e., with event detection, are highly recommendable in non-impulsive periods because 
of their high integration order. The drawback of event-driven schemes is that they cannot represent an infinite 
accumulation of impacts. On the other hand, timestepping schemes are applicable in impulsive periods but they 
are of integration order one in both impulsive and non-impulsive periods [21. Integration schemes known from 
computational mechanics ll28l[T7ll4^ l8l [T^I^ usually suffer from oscillations in the relative contact velocities 
at least if impulsive periods occur ETl . 

In |[3^ l40l l4T1l . classic timestepping schemes ll^ l25l |32l |33l are improved by splitting the non-impulsive 
and impulsive forces, such that the advantages of event-driven and timestepping schemes are conserved and 
the disadvantages are avoided. The approach is theoretically based on time-discontinuous Galerkin methods. 
Hence, we assume jumps for the velocity across discretization intervals, such that non-impulsive forces benefit 
from higher order trial functions and impulsive reactions yield local integration order one automatically. For 
the piecewise linear velocity trial functions discussed, an effective algorithm and framework are presented. In 
the case of finite element applications, it is shown however that the base integration scheme may get unstable 
without additional damping because of its half-explicit nature. That is why, we derive different timestepping 
methods based on the general framework of impulsive corrections and non-impulsive base integration schemes 
but leave the time-discontinuous Galerkin setting. We use base integration schemes like the generalized-a 
method IITSll . the Bathe-method [i61 and the ED-a method Q within mixed timestepping schemes EOl [H 
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for high-frequency damping instead. A comparison between acceleration and velocity level B |22l |T8l |42l 
for multi-contact cases is studied in BOl 1411 . We evaluate both the non-impulsive forces and the impulsive 
reactions on the advanced velocity level using the projection formulation for convex sets together with semi¬ 
smooth Newton methods ll4l [T4ll^l^ . 

Timestepping methods for structural dynamics are evaluated among other criteria based on two well understood 
and appreciated characteristics: unconditional stability and high frequency dissipation 10111241 [T5ll45l [8l. The 
first criterion ensures that the stability of the method does not depend on the time-step size. The other require¬ 
ment concerns the algorithmic dissipation and damping properties for non-physical high-frequency modes. In 
structural dynamics, large scale movements and lower modes are frequently of particular interest whereas high 
frequencies are not considered or only get into the system as spurious oscillations due to the discretization. 
However, the high-frequency modes affect the numerical and convergence properties of the system especially 
in the nonlinear regime. Hence, it is desirable to reduce these frequencies by numerical damping. The intensity 
of the damping should be controllable and should affect the lower modes as less as possible. In order to meet 
the above mentioned criteria, a variety of integration methods has been developed and optimized with regard 
to the desired properties. The HHT-formalism Il2^ and the generalized-a-method ifTSl are two well known 
methods based on the Newmark-scheme [[31] ■ Having a controlling parameter for high frequency damping and 
second-order accuracy for the generalized coordinates and the Lagrange multipliers mark the generalized-a 
method as a good choice for numerical integration [*5i|. An extended version of the generalized-a method for 
nonsmooth problems is presented in |[T3l [T2l . However, the global accuracy remains order one in the presence 
of constraint forces. An undesirable property is that the damping parameters have to be selected and acceptable 
values depend on the characteristics of the problem being solved. The damping parameter in Newmark-based 
integrators plays a very important role as the high frequency oscillations may effect the results and even con¬ 
vergence of the finite element solver |[T0ll40ll4T]| . Recently, new attempts have been made to overcome such 
problems, among those, we focus on two schemes: the Bathe-method and the ED-a method, which both use 
extra information in each time interval. The Bathe-method combines the use of the trapezoidal rule and the 
Euler backward method |I6|. It has been shown that the additional midpoint information can result in a more 
robust simulation in the nonlinear analysis and can avoid nonphysical large contact forces. The Bathe-method 
is also effective in the linear analysis [7]. The ED-a (Energy Decaying) scheme is proved to be one of the 
most effective methods in solving stiff nonlinear finite element problems. The ED-a method is based on a 
variational interpretation of Runge-Kutta methods and time-discontinuous Galerkin (TDG) methods |’9i It has 
been shown that the ED-a scheme benefits from unconditional stability in the nonlinear regime and also damps 
out the unresolved and spurious high frequencies using tunable parameters lITOl . 

Among different examples in the area of nonlinear structural dynamics, flexible mulfibody systems IMHH are 
of interest for the present paper. Sudden impacts and persistent contacts in addition to the nonlinear behavior of 
the system due to large displacements and deformation increase the difficulty to solve the system of equations 
properly and increase the level of sensitivity regarding stability 14011411 . As an example, we consider a flexible 
slider-crank mechanism with clearance and dry friction. The numerical damping capabilities are explored by 
means of an elastic connecting rod, while the other parts of the system are assumed to be rigid bodies. The 
derivation and evaluation of the system matrices and vectors is oriented towards Il44ll by applying the floating 
frame approach. Eor the elastic connecting rod, we consider a beam type element with three degrees of freedom 
at each node. Hence in Sect.|^ we explain the equations of motion of nonsmooth mechanical systems and mixed 
timestepping schemes, which benefit from high-frequency damping. Their damping behavior is compared for 
a simple linear system with large stiffness and bilateral constraint. In Sect. the slider-crank mechanism is 
introduced. Sect. carries out the validation of the algorithm by comparing the results with the rigid body 
studies in l|2T|. The convergence order is studied by comparing the results with a SimpacliQ simulation. The 
damping properties of different time integration schemes are discussed in detail based on the amplitudes of 
the frequency spectrum and percentage of period elongations. In addition, a modal analysis is undertaken. 
Sect. 1^ summarizes mixed timestepping schemes with high-frequency damping based on powerful and stable 
implicit time integration algorithms from computational mechanics such as the generalized-a method, the 
Bathe-method, and the ED-a method for complex nonlinear systems in structural dynamics. This paper is 
based on the student work lIMIl . 

^ http ://www. simpack.com/ 



2 Equations of motion and time-discretization 


In this section, we first derive the general form of the dynamic equations of motion for flexible multibody 
systems subject to unilateral contacts and friction (Fig. [T]l. An example is the slider-crank mechanism for a 



predefined gap between slider and the bordering wall. We discuss how such mechanical system can be modeled 
and how the time evolution of such systems can be obtained by numerical integration. Therefore, we introduce 
and discuss the properties of the generalized-a method, the Bathe-method, and the ED-a method for solving 
the equations of motion implicitly in time. For each method, we solve for impulsive forces A and non-impulsive 
forces A based on the idea of separation of impact and contact forces ll^ |40l |4T1 . 

2.1 Nonsmooth approach 

The nonsmooth approach is a method for modeling mechanical systems with unilateral contacts and friction 
using set-valued force laws. We consider the slider in Fig. [T] which slides, sticks or impacts on the cylinder 
wall. For each comer of the slider, normal and tangential gj gaps as well as normal and tangential Xj 
reaction forces have to be evaluated according to Fig.[^ 



Figure 2: Normal gap as well as normal and tangential contact force. 


The associated set-valued friction law of Coulomb-type is shown on the right-hand side of Fig.[^ Regarding the 
sliding case, the friction force is given explicitly. Regarding the sticking case, the friction force is set-valued and 
determined according to an additional algebraic constraint. The nonsmooth approach changes the underlying 










closed contact 


> 0,gN = 0 

-► 

gN 

open contact 

= 0,gN > 0 



Figure 3: Contact laws. 


mathematical structure if required and leads to a proper description of the mechanical systems (left-hand side 
of Fig. 1^. As a consequence of the changing mathematical structure, impacts can occur and the time evolution 
of the positions and the velocities cannot be assumed to be continuously differentiable anymore. Therefore, 
additional impact equations and impact laws have to be defined, which are discussed in the following sections. 

2.1.1 Timestepping integrators 

Timestepping schemes are based on a time-discretization of the system dynamics including the contact condi¬ 
tions in normal and tangential direction. The whole set of discretized equations and constraints is used to com¬ 
pute the next state of the motion. In contrast to event-driven schemes, these methods need no event-detection 
and are very robust in application. Moreover in the following section, a time discretization is used where the 
constraints are satisfied on velocity level. The accuracy of classic timestepping schemes is low. We propose 
higher order integration methods with high-frequency damping to process non-impulsive periods of the motion 
with larger time-step sizes and increased integration order. 

2.1.2 Normal contact law 

The contact law for normal constraints is depicted on the left-hand side of Fig. The normal contact force 
Xn vanishes if the bodies are separated > 0) and can only be positive if the bodies are in contact (gyy = 0). 
These constraints together form a complementarity condition: 

0<g^TAw>0. (1) 


2.1.3 Coulomb’s friction law 

Coulomb’s friction law is used to model friction forces. It states that the sliding friction force is proportional 
to the normal force of a contact. The amount of the static friction force is less than or equal to the maximum 
static friction force, which is also proportional to the normal force of a contact. The sliding friction force has 
the opposite direction of the relative velocity of the frictional contact. For closed contacts, the overall contact 
law in tangential direction is 


\\Xt\\<IJ-Xn for gr=0Ag^<0, 

for gj <0 


Xt = — 


11^7 


( 2 ) 


The coefficient of friction /r, in general, is a function of the relative velocity. We assume it to be a constant 
value. 









2.1.4 Projection formulation 


The projection formulation transforms the conditions shown in Fig. in equivalent equations using convex 
analysis. The proximal point to x in a convex set C is defined by the projection lITTII 

projc : IR —> IR , X proj(-(x) = argmin ||x —x*|| . (3) 

a:*6C 

With this definition, we express the normal contact law on velocity level as 

^N, - projiR+ (Aw, - rgyv,) = 0 (4) 

for all contacts belonging to the index set of closed constraints 

J^I = G J^o : gN, < 0} ; (5) 

where J^o contains all constraints. The arbitrary auxiliary parameter r > 0 represents the slope of the regular¬ 
izing function; it may be used for stabilizing the solution process lf39l . In the same manner, we formulate the 
Coulomb friction law Q as 

ProjcrfA^,) - rgr) = 0 (6) 

for all contacts in J’x , where the corresponding convex set is given by 

Ct : : y^Cr(y) = (x G IR^ | ||x|| < l^\y\} . (7) 


2.1.5 Newton’s law of impact 


For countable time instances tj, the evolution of the slider-crank mechanism might get impulsive. For some 
component k* of the gap function gWi* = 0 but gWj.* (^(0) > 0 for — 5 < t < tj and an appropriate 

5 > 0. This possibly leads to jumps in the velocity variables. Their derivatives do not exist anymore in the 
classical sense ll^l40ll4Tll . We have to define fhe leff-hand and righf-hand limits: 


Snj 

:=limg^{t) , 

tttj 

St, 

:=\imgT{t) , 
tVj 

(8) 

Snj 

:=limgw(0 , 

titj 

St, 

:=limgr(f) . 

t\ij 

(9) 


Then, the Lagrange multipliers describe the finite impulsive interaction in the sense of distributions: 

Awj = lim / Awdt, Ajj = bm / Aj-dt. 


5;o Jtj~s 


5^0 Jtj-S 


( 10 ) 


Newton’s impact law describes the elasticity of the collision by considering the local velocities before {gj ) and 
after (gt) the impact: 

0 < + SNgNi -L Aw, > 0 , 


I Ar, 11 < IJ-Anj for g^. -y e j.= 0 


St, 


Atj = - —I^Anj for g^. + ergr. 7^ 0 

St. 


( 11 ) 

( 12 ) 

(13) 


where Cw and Sj are the coefficients of restitution in normal and tangential direction, respectively. Therefore 
for the normal impact equations, we get 


Aw,- -projiR+ (Aw, - r{gjf^+eNgNj)) = 0 • 
For the tangential impact equations, we get 

Ar, -projc^(A„p (Ay. - r{g+ + ergy.)) = 0 . 


(14) 


(15) 







2.2 Equation of motion for a multibody system 

By determining the kinetic energy of the deformable bodies, the virtual work of the internal and external forces 
and the kinematic constraints that describe mechanical joints as well as specified bilateral or unilateral con¬ 
straints, one can use Lagrange’s equation to write the system equations of motion for multibody systems: 

Mq + Cq+Kq = h + WN)^N + yVT)iT , (16) 

where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, Wn and ILr are constraint 
matrices in normal and tangential direction, respectively, and Ay are normal and tangential contact forces, 
respectively, and h combines the effect of external forces and the quadratic velocity vector. If an impact occurs, 
the impact equations 


Mi 


vX — V ■ 
; j 


= + Ik Tj^Tj 


(17) 


Newton’s impact law ( [TT] ) with restitution coefficient G [0,1] and Newton’s impact law ( [T3| ) with Et G [0,1] 
have to be solved. 


2.3 Summary of computational algorithm 

The idea is to combine both the non-impulsive motion and the impulsive motion within one consistent integra¬ 
tion scheme. First, this integration scheme has to model impacts and velocity jumps automatically if necessary. 
Second, it has to switch to effective higher order integration with all kinds of nice benefits of sophisticated 
integration schemes for differential algebraic equations. We introduce a framework, which is derived from a 
time-discontinuous Galerkin setting as proposed in 1^1411 . It is more abstract and includes integration schemes 
for differential algebraic equations as base integration schemes on velocity level for one time-step from f, to 
?,+i. With this propagation at the end of the time-step, we can check by the same activity rules, i.e., index set 
calculations on velocity level with ^i, if new contacts have been closed: 

3k* : {qi) > 0 A (^,-+i) < 0 . (18) 

In this case, we just correct the solution for the velocity variables and calculate v^j with the impulsive forces 
Ayv,^i andAz^^j. If 


tk* : gN,, {qd > 0 A g^^, (^;+i) < 0 , (19) 

we just set The overall algorithm can be summarized as shown in Fig.|^ 

As an example. Fig. shows the slider at different time-steps before and after an impact. The dotted line 
represents the normal gap velocity before an impact. Going from to there is no new active contact 
point. Hence, we just have to solve for non-impulsive forces on velocity level to avoid further penetration. We 
notice the usual drift-off effect (gA^^ (^z+e) — (^) imposing the contact forces A a? on velocity level which 
only constrains g^.^^ > 0. 

We proceed with explaining the generalized-a method as a base integration scheme. Then, we focus on how 
to calculate impulsive corrections. The Bathe-method as well as the ED-a method are further base integration 
schemes, which can be used instead of the generalized-a method. We explain them and compare the properties 
of all three base integration schemes at the end of the following section. 




specify characteristics, end time T, 
time-step size At 


t = 0 



2.4 Generalized-a method 

According to the generalized-a method for flexible multibody systems can be summarized as follows: 


+C;+iVj-^j + — 

, ( 20 ) 
(1 + CCffi^i — (1 “b CCf^i j (21) 

‘ii+i =^/+Atvb +At^ [(0.5 — ^)A,-|-^A,-+i] , (22) 

+ A? [(1 — 7)A,'+ 7A,'+i] , (23) 

% = q{^) , = V(0) , (24) 

do =Mo^ (/jo-i^o^o“<^ovo) , Ao = ao, (25) 


where is the vector of generalized coordinates, v,- is the vector of generalized velocities, a,- is the vector 
of generalized accelerations. A, is the vector of acceleration-like auxiliary variables, defined by the recurrence 
relation (|2T]), and At is the time-step size. 
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Figure 5: Transition from impulsive to non-impulsive reactions. 


2.4.1 General characteristics 

The displacement and velocity update ( |2^ and ( |2^ are identical to those of the Newmark algorithm. The 
structure of these update equations is obtained using Taylor series expansion about f,. The crucial task is to 
determine the relationship between the algorithmic parameters, am, af, 7 , and j3. With appropriate expressions 
for 7 and p, and if a„, = 0, the algorithm reduces to the HHT-a method. The generalized-a method, with 
parametric values given in ( | 2 ^ , is unconditionally stable for linear problems, second order accurate possessing 
an optimal combination of high-frequency and low-frequency dissipation: 


OJm - 


«/ = 


2poo - 1 
Poo 1 

Poo 

Poo 1 
1 


(26) 


7 = 2 “ ’ 

1 2 

^ = — (1 — tt,j, -|- af) , 


where Poo G [0,1] is the spectral radius of the amplification matrix at the high frequency limit. The stability 
region is indicated by the shaded area in Fig.[^(|^. In order to have a better insight on the effect of the different 
parameters on the eigenvalues of the amplification matrix, we plot the spectral radius and the relative period 
error with respect to At/T, i.e., the time-step divided by the period T = 27ij(0. First, we consider a special 
case of the generalized-a method by setting a^ = 0 (HHT method). From Fig. we are allowed to choose 
af between 0 and 0.5 (Fig. |^. We see a cusp after we pass the point E. Using optimum values for the 
generalized-a method as introduced in (2^\, we see the behavior of the method for different poo in Fig.j^from 
the no-dissipation case (Poc=l) to the so-called asymptotic annihilation case (poo= 0 ), or moving along the red 
dotted line from point D to point B. In order to observe the properties of different regions in Fig.[^ we plot the 
spectral radius for point A and point C (Fig. |^. If we select the a„, and ay values away from the dotted line 

(||A“|| = A “2 ), we may expect the cusp. Figure 9 shows how we can modify the a,„ and a/ values for point 
A and point C in order to have the same p<x, but witn smooth transition instead. Further properties can be found 
inlAl 
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Figure 6: Generalized-a stability region in a,n-(Xf space, different test case A(—1,1/2), B(—1,0), 
C(-l, -1/2), D(l/2,1/2), £(0,1/3) O- 

2.4.2 Calculation of contact forces on velocity level 

Using ([T]) and Q together with ( [T^ , we calculate contact forces and We interpret ([T]l and Q on 

velocity level |[4^ |4T|| : 


- 


- 


if gNi > 0 


,proJiR+ 

'O 


■^M+i '' gNi- 


/+1 


else 




^T.^i-rgn 


/+1 


if gNi > 0 
else 


(27) 

(28) 


using the projection-formulation row-by-row and a predictor for closed contacts gjq. < 0. Thereby, the local 
velocities satisfy 


gNi^i = > 




(29) 


and depend on and XTi^^ because of the implicit representation. Hence, we calculate equivalent forces 
and such that the given local velocities are projected into their respective admissible space. Using 
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Figure 7: Properties of the HHT-method for different a = af, am = 0. 


( |2T] ) and (^i, we substitute: 




Sn,+, = ^Ni+i +Ary _ J ] , 

1 CXm 


_ J {Gtn,+ i>^N,+x +GTi^^XT,^^ } , 

i IXm 


where 


^^,.1 = + At (1 - 7) 


1 — ttf ^-1 ^ af am 


'+' V 1 - a, 


1 OJm' 


Ft,^. = vt + At (1 - 7) Wl^^M + AtyWl 


1 OJ/71 

an 


^ R \ n 

Ali+iKi+l + , _ ] ) 

i IXm; 1 IXti- 


+‘ V 1 - a, 


-A, 


and 

, Gmt.,, = > 

are called Delassus matrices with the effective mass matrix 

Mi+i =Mi+i +Atiy\^Ci+,+At}^\^Ki+i 

1 am 1 am 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


and the effective right-hand side 


an 


«/ 


Ri+1 = - Q+i V; - [ At,' (1 - 7) - At,'7 ) C,'+iA,' - At,'7-——Q+ia,- - ^z+i^,- - At,'/:,'+i v4 (37) 

\ t OCm / t On 


- Atf{0.5-p)-Atfp 


an 


1 


Ki+iAi-Atfp 


«/ 


I-a. 


-Ki+iai 


(38) 


We focus on active contacts and transform ( [27| ) and ( [28| ) formally using row-by-row interpretation: 

A-7 


"7]+i,- projcrr;+i7;-+,' 


(40) 
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Figure 8: Properties of the generalized-a method for different poo. 


Thereby, the index set of closed constraints at r, is given by 

J^I = {k£J^0 : gt^^{qi) < 0 } . 


(41) 


In the multi-contact case, active contacts might be depending. Hence, we cannot use a nonsmooth Newton 
method but we solve (|^i, ( |40| ) by a nonsmooth variant of the Gauss-Newton method, i.e., we choose an 
approximate root (Aat j?,) of the function 


/ : IRI‘^‘1 xIR- 


IR '^' xIR 


2 |j^l 




X 


- proj]R+ 


- fgN,A 


-projcr(Ajv„,j) ^T,yi -rgj j;^{XN,J!i,XT,,^^) 


(42) 


with the Moore-Penrose pseudoinverse operator pinv. Using the Gauss-Newton algorithm, we calculate un¬ 
known contact forces with the following iterative algorithm: 

= ( 0 , 0 ) , / = 

while ||/|| > tol 

; \ ®N{-rAty\:^GNT,j?i) 

T, ) “ PIWV ( V/(A jT, , A jTj ) )/ 

/ = / {Xn,J!i,Xt,.^i) 

end 

The occurring Heaviside functions 




0^ : IrI'^'I X IR^I^il ^ diagl-^‘l’l^‘1 , 0^,, (A^,^,, Ar,^,) 

0r : IrI^‘1 xIR2|^iI ^ diagl=^‘l’l-^‘l , 0r,, {Xn,^„Xt,^,) 



if Aiv, - rgf^^ < 0 

1 

(44) 

else 


if 

K - rgT, 

> p AwJ 

(45) 

else 



are interpreted row-by-row and we use r = 0.1 without adaptation to improve convergence of the numerical 
scheme ||39l. 
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Figure 9: Smooth transition of cases A and C by modifying af and 


2.5 Assignment of impulses 


If on the other hand, there is at least one inactive gap function (gyv*. {‘li) > 0) at U that becomes active {gN,^ i^i+i) < 
0), an impact occurs in a rigidly connected component of its source. Hence, we update the velocity at ti+i using 

O: 

= Vi +^m^M+i^W/+i + T.+Ar.+i , (46) 


where ( ) and ( denote the parameter at tj+i before and after the impact, respectively. For the computation 
of v^j, Aa?. and Aj^, we write ( [TT] ) and ( pA] ) on velocity level row-by-row: 


Aw,,, = < 


if > 0 


proJiR+ 


Aw,+,+ewgyv,^,) 


SNi 


1+1 


else 


At^+i - < 


P™JCAA„,^ 


^T.^i-r{gT,+,+eTgT^J 




if§M+, >0 


else 


Equation ( [TT] ) can be used to eliminate and g^ , row-by-row resulting in 




A7;._n,Jf/+i-Pl‘ojcr(A ,+i) 


(47) 


(48) 


(49) 

(50) 


Again we look for the roots of 

/ : ^IrIAI xIR^I'^‘1 

(Aw.j^hAt’^j^j) /{AnAt, j^i) = 


( Ayv,- proj]j^+ [Ayv,Jf, - rgyv jr, {An^j^i , Aj^jTi )] 

-projcr(A^.,|) [Ar,.^i -rgT,yi{-^N,A,-^T,yi)] 


( 51 ) 
































with a nonsmooth Gauss-Newton method as in 


( |4^ . The derivative of / at (A^v.j^i jr,) is given by 


// — 0 A?(/ — &Ni—rGNT,j^i) ^ 
&Ti—rGTN,j^i) I — ®t{I — rGT^,^i)) 

The occurring Heaviside functions 



@N 


@T 


IrI-^.I xIR 2|A| ^diagl'^‘l’l=^‘l (A^v,^,, A^,^,) = 

xIR^I^iUdiagl-^'I’l-^'l @T,, {An,At,a) = i 


if An, - rgN, < 0 
else ’ 

if ||A7i-rgrJ| > ft |Aw* 
else 


( 52 ) 


(53) 

(54) 


are interpreted row-by-row. We use r = 0.1 A/ without adaptation to improve convergence of the numerical 
scheme. In nonlinear dynamics, we have to apply an iterative algorithm to calculate unknowns at f,+i based on 
updating mass, stiffness, damping, gap functions and external forces. In case of no impact applying the fixed 
point iteration method, we obtain the unknowns at using the generalized-a algorithm for each time step: 


Generalized-a time integration scheme 


Set yt = 0 

<?°+i = <?/> ■^’?+i = A-, a^i+i = ai, Af+i = A; 

Ml, = Mu hi, = hi, = Wn., = W% Kl, = Ki, cl, = Q 

while true 

calculate with Gauss-Newton algorithm 


(7^+1 — 

“;+i ~ 


M,. 


(+1 


" it 


A^+l — ^ 1 I 

'+1 “ 1 ^'+1 + 


jk 

CCf OLn 


;+l +'^Ni+iAN,+i 




-a, — 


-A,- 


l-a„, ^ l-Um l-a,r. 

V •+/ = Vi + Ati{l-Y) Ai + AtiyAll 

‘?f+i = + AtiVi + Atf (0.5 — [5)Ai + Atf pA^ 


'■+1 


< tol = vH break 


Mil =M [kill 

dfi'=c(9‘;;,v;+'), K‘*; = k I+ia;) 
"'ft', = »'» (9m ). wSl = 'fr l\ 

k = k+ 1 


A+i 


,hll=h 


^!+l A /+1 


end 


(55) 


2.6 Bathe-method 

The Bathe-method fT] is an effective implicit time integration scheme, which has been proposed for the finite 
element solution of nonlinear problems in structural dynamics. Various important attributes have been demon¬ 
strated. In particular, it has been shown that the scheme remains stable without the use of adjustable parameters. 
For this method, the complete time step At is subdivided into two equal sub-steps. For the first sub-step the 
trapezoidal rule is used and for the second sub-step the 3-point Euler backward method is employed, as it is 
described in ( |58| ) to ( [^ . 











2.6.1 General characteristics 


According to [7], we have to consider the dynamic equilibrium for time t + At and t +At/2 which is indicated 
using index i + 1 and index / + 1/2: 


l -\-2 '^2 ^^2 '^^2 *^^2 ^^2 i+j i+j i+j i+j 
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(56) 

(57) 

(58) 

(59) 

(60) 

(61) 

(62) 


In Fig.[^for poo = 0, we compare the Bathe-method and the generalized-a method. We notice that the Bathe- 
method preserves the low-frequency oscillation better than the similar generalized-a method, and of course 
both avoid high-frequency vibrations. Total annihilation is important to avoid high frequency noises generated 
by the time integration algorithm. The Bathe-method shows a better behavior for the period error in comparison 
to the generalized-a method and the classic Newmark method. Further properties can be found in|^ 




At/T 


Figure 10: Comparison of second-order methods: Newmark method, generalized-a method and Bathe-method. 


2.6.2 Calculation of contact forces on velocity level 

Using the Bathe-method (|58|)-([^, we have to modify the implicit representation of gap velocities in ( [29l ). We 
have to solve for the unknowns first at t + At/2 and then at t -I- At: 
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After calculation of the unknowns at t,+i /2 for the second half of the interval, we have: 
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The matrices Ki, K 2 , Ri and R 2 are defined in (|75)-(|76|). Once we have calculated the velocities v^^j, the 
computation of the impulsive forces is given in Sect. 2.5 The generalized-a method and the Bathe-method 
use the same procedure to update the calculated velocities after the impact. In case of no impacts, applying the 






fixed point iteration method for each time step of the Bathe-method, we get: 


Bathe time integration scheme: first half step 
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Bathe time integration scheme: second half step 
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(76) 



2.7 ED-a method 


In energy decaying schemes, we develop robust algorithms for integrating stiff nonlinear finite element prob¬ 
lems in time. The basic motivation behind these schemes is that classical algorithms that are unconditionally 
stable and high frequency dissipate in the linear regime, loose their properties in the nonlinear regime ifT^ . We 
typically interpret the stages qj, vj, aj as field variables associated with the time In this sense, the unknown 
fields are allowed fo creafe a jump discontinuify af fhe beginning of fhe time sfep thaf is responsible for the high 
frequency damping behavior of the scheme f9j. 


MjQj + CjVj+Kjqj = hj + WTj^Tj , 
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^,+1 =^, +Y , 

flo = Mq ^ {ho — Koqo —C). 


(77) 

(78) 

(79) 

(80) 
(81) 
(82) 
(83) 


In this work for practical implementation of the scheme, the displacements qj and q^^^ are eliminated, leaving a 
velocity-based iteration scheme in the 2 x ndof unknowns vj and v,+i. However, the overall procedure is more 
expensive than other one-stage schemes like the generalized-a method, since the matrices are twice as large. 


Note that for Uar = 0 or a = 0, we recover a conserving scheme. The parameter War does not control the 
asymptotic value of the spectral radius but only controls the cut-off frequency of the scheme and so relative 
period errors (Fig. 111. The minimum period elongation is obtained for War = 1/6. The method is second order 
accurate for arbitrary War > 0 and arbitrary ordinary differential equations; third order accuracy is obtained for 
the scalar linear model problem in Sect. 2.8 and the special value War = 1/6. The parameter a is responsible 
for the asymptotic value of the spectral radius, in fact ((84l) is an optimal choice ifTOll (Fig. 12 1 : 
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In case of no impact using (f/T]) to ([82l) we write: 
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Figure 11: Comparison of methods: ED- 



MjT 

method with Poo = 0 and Bathe-method. 


ED-a time integration scheme 
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The matrices Me, Rc, WWj^ and vectors Vc, Ac are defined as 
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Figure 12: Comparison of methods: ED-a method with Uar = 1/6 and generalized-a method 
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Figure 13: ED-a method in time. 
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2.8 Comparison of base integration schemes 

Our objective in this section is to present the solution of a simple linear system as a model problem to represent 
the stiff and flexible parts. We compare the presented base integration schemes and discuss their properties 
within the overall framework. We also compare the results for bilateral contact forces using the acceleration 


























level and velocity level approach to satisfy the constraints. Let us consider the solution of the 3 degree of 
freedom mass-spring system shown in Fig.[^for which the governing equations are 


/mi 0 0\ /qi\ /ki +k2 

0 m 2 0 ^2 + “^2 

\ 0 0 ma/ \^3/ \ 0 


-k2 0 \ /mY\ 

k2 + k3 -kj, ^2 = my 

-k-i ^3 + ^ 0 / \? 3 / \my) 


(91) 


We use the parameters /ci = 1 N/m, ^2 = 1N/m, kj, = 1 N/m, ^0 = 10^ N/m, mi = 1 kg, m 2 = 1 kg, m 3 = 1 kg, 
7 = 9.81 m/s^. In the left part of Fig. 14 a stiff spring is used to represent, for example, an almost rigid 


connections Q, while the other springs represent the flexible parts of the structural model. 




Figure 14: Model problem of a three degree of freedom mass-spring system, (a): model with stiff spring, (b): 
model with bilateral constraint. 

The highest frequency in case (a) is /max = 503.3 Hz, which is due to the stiff spring with stiffness ko- The exact 
trajectory for mass number 3 is plotted using At = 10^^ s to represent it more properly, and the time-step size 
for comparison of different base integration schemes is chosen as At = 10^^ s. Figure [T^ shows the response 
of the first two degrees of freedom which are connected with the soft springs. As we observe, the response 
from different base integration schemes follows the reference solution coming from modal analysis. The only 
important difference is the phase shift due to numerical integration which is minimum for the ED-a method. 
The response of mass number 3 which is connected to the rigid wall with a stiff spring is shown in Fig. [T^ 
The generalized-a method with poo = 0.8 is able to represent the high frequency vibration with some phase 
shift, however the other methods damp out this high frequency mode in a few time steps. Using At = 10^^ s, 
i.e., a sampling frequency fs = 1000Hz, we have to make sure that all frequencies higher than ^ = 500Hz are 
damped out, in order to avoid chattering in the response. 

Figure [T7] shows the difference between simulation results if we try to satisfy the constraint on acceleration and 
velocity level. We conclude that using the acceleration level results in a constant residual velocity for mass 
3, and of course a linear drift-off effect and violation of the constraint. The calculation of contact forces is 
straightforward and follows the reference solution (Fig. [T8] ). Using the velocity level, the contact force needs 
to be calculated with an additional integration. This calculation results in a constant position error for this 
model problem. Using the velocity level results in an oscillatory behavior for the acceleration of mass 3 which 
means that we expect the same oscillatory behavior for the calculation of contact forces (Fig. [T^. Inserting 
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Figure 15: Displacement, velocity and acceleration of masses 1 and 2 for various integration methods. 


additional damping for the numerical integration results in better approximations of contact forces and avoids 
those artificial high frequencies. The results using the Bathe-method or the ED-a method are similar to the 
results using the generalized-a method with poo = 0.0. It is worth to mention that the highest frequency in the 
system is /max = 0.8717 Hz, i.e., ^ = 0.0436, which means that the oscillatory behavior of the contact force in 
case of no damping comes from the structure of the generalized-a method and has nothing to do with the poor 
representation of high frequencies in the system (model (a)). We conclude that the unsymmetrical structure of 
the Bathe-method or the ED-a method helps to improve the calculation of contact forces in case of model (b). 

A similar observation has already been formulated in HOl IdTl for the application of half-explicit timestepping 
schemes on velocity level on examples, which discuss an impacting elastic bar or a rubbing rotor. 


3 Application to a flexible multibody system 


We consider the slider-crank mechanism shown in Eig. 19 where h is the length of the crank and I 2 the initial 
length of the connecting rod for the undeformed state. The inertia of rigid crank, flexible connecting rod and 
rigid slider consist of translational masses, i.e., mi, m 2 and m 3 , as well as of rotational inertia values, i.e., Ji, J 2 
and 73 for the undeformed state. Eor the flexible connecting rod, we consider p and E as density and Young’s 
modulus, respectively. The cross-sectional area is given by A = HD, i.e., the product of height H = m 2 /{pl 2 D) 


and thickness D = y \ 2 J 2 lm 2 —of the rod, and the second moment of area is given by 7 = . The 

system is subject to gravitation 7 Ell¬ 


in order to describe the flexible system, the local coordinate system of the connecting rod is located tangentially 
in the joint between the crank and the connecting rod. Thus, the basis for a floating frame of reference formu¬ 
lation is accomplished. Such a description is characterized by a separation of the coordinates of an elastic body 
into reference and elastic coordinates. The reference coordinates delineate the rigid body movement and consist 
of the translational coordinates describing the absolute position of the local coordinate system and the rotational 
coordinates describing the orientation by angles. The elastic coordinates capture the flexible movement. Crank 
and slider are described by minimal coordinates. Eor the evaluation of the equations of motion, we consider the 
inertia coupling between the different sets of coordinates. The derivation of the system matrices for the given 
slider-crank mechanism is studied inICl 




















Figure 16: Displacement of mass 3 for various integration methods. 


4 Results 

The application of the generalized-a method, the Bathe-method and the ED-a method as base integration 
schemes of the overall framework to nonsmooth problems with unilateral contacts with friction is discussed in 
the following section. First, we show the spatial convergence of the schemes and validate the results comparing 
to the rigid case. Finally, we discuss the time integration schemes concerning different aspects. 


4.1 Validation for a rigid slider-crank mechanism 

Figure shows the convergence for the position of the slider mass center when we increase the number of 
elements and run the simulation with At = 10^^ s using the Bathe-method. Specific characteristics are given 
in Table The convergence results are comparable for all base integration schemes. For the validation of the 
results, we compare an almost rigid system {E = lO^^N/m^) with a rigid multibody system with 3 degrees of 
freedom. Figurej^shows the results for the position of the slider mass center compared to the same simulation 
with rigid bodies 1211. The generalized-a method and the Bathe-method are both second order accurate. In 
order to show this behavior in a simulation, we consider the specific case of a bilaferal confacf by seffing fhe gap 
c = Om. Figure 1^ shows fhe convergence when we decrease At. The Slope of fhe line m shows fhe accuracy 
of fhe calculations in logarifhmic scale. The reference solufion is calculated wifh Simpaclj^ 


4.2 Comparison between generalized-a and Bathe-method 


In order fo have a heller insighl for closed gap silualions and fhe calculalion of confacf forces, some inilial 
properlies of fhe slider-crank mechanism are changed according lo Table whereas fhe olher characferisfics 
are sel according lo Table[^(Fig. [2^ . The lime-step size is At = 10^^ s. In Fig. 24 we see the comparison of 
the schemes for the angular displacements 0 i and 62 , i-e., the inclinations of crank and connecting rod, (63 = 0 ). 
Both algorithms with different damping values behave almost the same, before and after the impact. Figures [25] 
and [^ show the comparison for angular velocity and torque at the beam root. If we do not consider enough 
damping for the time integration schemes, high frequency oscillations corrupt the system response, especially in 
the velocity and stress fields. In fhis particular case, fhe resulls for fhe generalized-a melhod wifh = 0.8,0.5 
gel unslable soon afler t = 0.05 s. This algorifhm Iransfers energy from fhe higher (arlificial) fo fhe lower 
meaningful modes. Figure [^ shows fhe comparison of fhe confacf force using fhe velocily level approach. 


^ http://www.simpack.com/ 
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Figure 17: Displacement, velocity and acceleration for mass 1, 2 and 3. 


The Bathe-method is able to represent the vibration as the contact is closed, whereas the generalized-a method 
results in an additional vibration for the velocity which causes the contact condition to go on and off repeatedly 
in time. 

4.3 Comparison between ED-a and Bathe-method 

Using the data in Tablej^with time-step size At = 10^"*s, we run the same simulations up to t = 0.5s to show 
the stability and robustness of the ED-a method and the Bathe-method. Figure shows the comparison of 
angular velocities which match perfectly. The phase shift can be explain with the relative period diagram in 
Fig.[TT] The effect of the phase shift is also noticeable in the calculation of the normal contact force in Fig.[29| 
The general behavior and amplitude of the contact force is in good agreement using energy decaying methods. 
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Figure 18: Comparison of acceleration and velocity level for the calculation of contact forces. 
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Figure 19: Slider-crank mechanism. 














y-position 


Geometrical characteristics 

l\ = 0.1530m (length crank) 
h = 0.3060m (length rod) 
a = 0.0500 m (half-length slider) 
b = 0.0250 m (half-height slider) 
c = 0 . 0010 m (gap) 

Inertia properties 

m\ = 0.0380kg (mass crank) 
m 2 = 0.0380kg (mass rod) 
m 3 = 0.0760 kg (mass slider) 

Ji = 7.4- lO^^kgm^ (inertia crank) 

J 2 = 5.9• lO^^kgm^ (inertia rod) 

J 2 = 2.7 • lO^^kgm^ (inertia slider) 

Force elements 

7 = 9.81 m/s^ (gravitation) 

Contact parameters 

evi = £Ni = £Ni = EVt = 0.4 

for slider comers 

£Ti = £t2 = £ 7 , = £Ti = 0.0 

Ml = /^2 = Ms = M4 = 0.01 

Initial conditions 

die = 0.0 

02„=O.O 

03 „=o.o 

£04 = 150.0rad/s 

0 > 2 f, = —75.0rad/s 
a) 3 „ = O.Orad/s 

Material properties 

£ = 2-10“N/m^ 

of flexible rod 

p =7800kg/m3 


Table 1: Characteristics of the slider-crank mechanism with unilateral constraints and friction ll^ . 
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Figure 20: Mass center movement for different numbers of elements. 


Geometrical characteristics 

c = 0.0005 m 

Driven torque for crank 

T = 1 N/m 

Contact parameters 

£n, = £n2 = £n3 = £Ni = 0.1 

for slider comers 

Ml =M2=M3 = M4 = 0.1 

Initial conditions 

a)i„ = O.Orad/s 


a> 2 „ = O.Orad/s 


a) 3 |, = O.Orad/s 


Table 2: Modified characteristics of the slider-crank mechanism. 
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Figure 21: Mass center movement; comparison to a rigid body simulation. 
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Figure 22: Relative errors of position, velocity and contact force. 





Figure 23: Slider crank modified configuration. 
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Figure 24: Comparison between generalized-a method and Bathe-method for the angles di and 62 - 



Figure 25: Comparison between generalized-a method and Bathe-method for the angular velocities G)i and ft) 2 . 
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Figure 26: Comparison between generalized-a method and Bathe-method for the torque at the beam root. 
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Figure 27: Comparison between generalized-a method and Bathe-method for the active component of the 
normal contact force. 
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Figure 28: Comparison between ED-a method with p<x> = 0 and Bathe-method for the angular velocities. 




Figure 29: Comparison between ED-a method with p„ = 0 and Bathe-method for the normal contact force. 

































































































4.4 Comparison of computing time 


Based on the setting in Table [T] we analyze the relative central processing unit (CPU) time for the computation 
of a rigid slider-crank mechanism. Thereby, we compare the generalized-a method, the Bathe-method and the 
ED-a method also with Moreau’s midpoint rule, which is a classic timestepping scheme, and the half-explicit 
timestepping scheme (HETS) proposed in BOlldTIl . 

Eirst, we regard the relative CPU time per time-step, exemplary for At = 10^^ s in Table for the bilateral 
case. We compare the necessary time-step sizes At and their corresponding relative error with respect to the 


Moreau 

HETS 

generalized-a 

Bathe 

ED-a 

Rel. CPU time 1.0 

1.35 

4.15 

4.35 

5.00 


Table 3: Relative CPU time for At = 10 ^ s (bilateral). 


reference Simpack solution in Table Thereby, we compute the relative error examplary for the connecting 
rod inclination 62 for a set of considered time instances 


err = 


\Q2{tk) - Q2,a{tk)\ ___ 
\^2,a{tk)\ 


(92) 


Eor a given time-step size, the computational effort of the new methods is minimal using the generalized-a 
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04 
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Table 4: Comparison of the error e for different At (T = 0.05 s, bilateral). 


method as there is no additional midpoint calculation. Eor the Bathe-method and the ED-a method, we use 
one additional point in the time integration algorithm which explains the increase in the computing time. In 
case of the ED-a method, the additional point can be interpreted as a jump at the beginning of the interval and 
according to ( [85] ), we increase the unknowns which have to be solved simultaneously by a factor of 2. The 
Bathe-method selects the additional point in the middle of the interval. The unknowns are solved independently 
from the unknowns at ti+\. Smaller relative time-steps (y) result in smaller changes of the mass matrix, stiffness 
matrix and force vector and thus in less calculation time for the unknowns compared to the ED-a method as 
we need less iterations to find the unknowns at the middle of the interval compared to the end of the interval. 
The old methods are much faster per time-step because of low-order or explicit evaluations. However for a 
complete comparison, we have to consider also Table Eor the same error, e.g. 10^', the time-step size for the 
classic timestepping has to be chosen much smaller than for the remaining four schemes. The benefit of the new 
schemes in comparison to the half-explicit timestepping, which is also of second order, is the high-frequency 
damping. 


Eor the unilateral case, the solution of the ED-a method with = 0.0, At = 10 ^ s is chosen as the reference 
solution. Results are given in Tables [^ and They confirm fhe resulfs of fhe bilaferal case and are even 


Moreau 

HETS 

generalized-a 

Bathe 

ED-a 

Rel. CPU time 1.0 

1.15 

2.72 

3.15 

3.40 


Table 5: Relafive CPU lime for At = 10 ^ s (unilateral). 


advantageous for the new methods because of their increased stability. 
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to 


10 

-3 s 
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-6 s 

Moreau 

6.7- 


6.4- 

10-3 

6.4- 

10-4 

HETS 

1 .1- 

10-3 

2 .8- 

10-3 

2.3- 

10-6 

generalized-a 

2 .0- 

10-3 

1 .2- 

10-4 

1.4- 

10-3 

Bathe 

1.4- 

10-3 

7.8- 

10-3 

8.9- 

10-6 

ED-a 

9.4- 

10-4 

2 .8- 

10-3 

2.3- 

10-6 


Table 6: Comparison of error e for different At {T = 0.05 s, unilateral). 


4.5 Modal approach 


A major advantage of the floating frame of reference formulation is that the finite element nodal coordinates 
can be easily reduced using modal analysis techniques, based on a reduced set of eigenvectors ©■ 


The spatial convergence is plotted in Fig. 30 and Fig. 31 when we increase the numbers of mode shapes. 
Thereby, the simulation is based on Table the time-step size is At = 10^^ s, the boundary condition is con¬ 
sidered as clamped-free (tangential), and we use the Bathe-method for time integration. After an impact oc¬ 
curred, depending on the boundary condition and impact points, we need more number of modes to describe 
the vibration behavior than in the non-impulsive period (before the impact). As far as we deal with nonlinear 
problems, slightly different conditions before the impact may result in large changes in time after the impact, 
which explains the differences in the results. Considering spectral analysis for the Bathe-method (Fig. [T0|) af¬ 
ter ^ = y = 10^, there will be no effective mode in the response. We simply calculate a sufficient number 
of modes considering the high frequency dissipation of the Bathe-method. In this case, all frequencies up to 
fc = = lO’ 1/s (which covers the first 52 modes out of total 63 modes) can be effective in the final so¬ 

lution. Results are in good agreement between modal and full FEM simulation for the calculation of normal 
contact forces. 


xIO"* xIO’'^ 




Figure 30: Mass center portrait for modal approach and tangential boundary conditions. 


Figure|^shows the comparison between modal solutions for different boundary conditions. Different boundary 
conditions result in different mode shapes and frequencies. Assuming the same impact point and condition we 
get different vibration behavior, but the general pattern may be similar. 
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Figure 31: Normal contact force convergence for modal approach. 


5 Summary and Conclusion 

This work deals with the consistent and efficient integration of nonsmooth flexible multibody systems with 
impacts and dry friction. We present a timestepping scheme, which evolves from the idea of time-discontinuous 
Galerkin methods ll40l |4T]| . However, we abstract this origin and develop a framework which improves a non- 
impulsive trajectory by impulsive correction after each time-step if necessary. This correction is automatic and 
is evaluated on the same kinematic level as the piecewise non-impulsive trajectory, i.e., on velocity level. The 
resulting overall mixed timestepping scheme is consistent for impulses and benefits from higher order in non- 
impulsive periods and all advantages of the base integration schemes used to calculate the approximation per 
time-step. 

We present a nonsmooth adaptation of the generalized-a method, the Bathe method and the ED-a method. It is 
applied to a slider-crank mechanism with a flexible connecting rod, impacts and dry friction. The elastic behav¬ 
ior of the connecting rod is compared using the different base integration schemes and a modal approach. The 
results are validated with respect to ll^ and a rigid body simulation in Simpack. It is shown, that introducing 
enough damping for Newmark-type integrators like the generalized-a method results in more stable solutions 
for nonsmooth multibody systems especially on velocity and acceleration level in comparison to BOl l4ll . It 
seems that schemes like the Bathe-method and the ED-a method behave more robust in the nonlinear regime. 
These methods are more expensive per time-step but less steps can be used and the methods remain stable, 
even if the Newmark-type integrator fails for large deformations and long time duration dynamic response 
calculations. 

Eor the future, it is valuable to pursue deeper mathematical analysis of the overall framework with different 
base integration schemes to prove the characteristics observed numerically especially for nonlinear multibody 
simulation. 

A Generalized-a method 

Eor the analysis of the generalized-a method, it is advantageous to reduce the coupled equation of motion to a 
series of uncoupled single degree of freedom systems using modal analysis and using eigenvector orthogonality. 
The linear single degree of freedom system with angular frequency ft) is given by 

ij + co^q = 0 , 


(93) 



















































Figure 32: Comparison between FEM and modal solution for different boundary conditions. 


where the terms related to external damping and forces are set to zero to study accuracy and stability properties 
of the algorithm. The generalized-a method described in ([20l) to (|25]) can be written in the compact form 


Xi+i—AgaXi, / G {0, 1,..., — 1} , 


(94) 


where X, = (< 7 ,-, Afv/, At^a,) and Ag^ is the amplification matrix for the generalized-a method. With 

D = 1 - ttm + (1 - a/ ) , 

D. = (oAt , 
ft) = 


(95) 


Im 


it is 


Aga — — 


^ /l-am-OCfpQ^ 


D 


— 

-Q? 


1 a,jj (2 (1 ®m) 

I-am-{I- a/) ( 7 - j8)a2 (1 - 7 ) (1 - a,„) - yUm - (1 - a/) (| - j8) D? 
-{l-af)D? - (1 - a/) (2 -jS) - a„, 

(96) 


The accuracy of an algorithm can be determined using the difference equation in terms of the displacement 

qi+i -Aiqi+A2qi-i -A3^,_2 = 0 , (97) 


where Ai is the trace of Aga, A 2 is the sum of the principal minors of Aga and A 3 is the determinant of Aga- It 
can be shown IITSl that the algorithm is second order accurate for unconstrained mechanical systems if 


1 

7 = 2 ““«« + “/• 


(98) 


The parameter 7 is responsible for numerical dissipation. For 7 = 0.5, there is no numerical dissipation, whereas 
for values 7 > 0.5 the numerical dissipation increases. The stability and numerical behavior of an algorithm 
depends on the eigenvalues of the amplification matrix. The spectral radius p of an algorithm is defined by 


p =max(|Ai|, IA 2 I, IA 3 I) , (99) 

where A, is fhe ifh eigenvalue of Aga- An algorifhm is uncondifionally sfable for linear problems if p < 1 for all 
G [0,oo). The generalized-a mefhod is uncondifionally sfable provided 

1 oil. 

«„,<«/<-, j 8 > - -h - (a/ - am) ■ 


(100) 

























The spectral radius is a measure for numerical dissipation. A smaller spectral radius corresponds to greater 
numerical dissipation. Desirable dissipation properties have spectral radius close to unity in the low-frequency 
domain; the value smoothly decreases as Q. increases. Typically, in the low-frequency domain, IA 3 I < |Ai, 2 |- To 
preserve the smoothness as Q. increases, IA 3 I < IA 1 . 2 I for all Q. G [0,oo). Violation of this condition will result in 
a cusp in the spectral radius plot where p increases as Q. increases (point A in Fig.[^. Calculating lim^i^coAga 
in (96 1 , we get the following eigenvalues of the amplification matrix in the high-frequency domain: 


A3“ = ^, 

^ ttf-l 


( 101 ) 


where j = \/—1. High-frequency dissipation is maximized if the principal roots (^“ 2 ) become real, i.e., 
3 (^“ 2 ) = 0. It can be shown from (|101|) that for the generalized-a method, this condition is satisfied if 


j8 = -{l-a^ + afY 


( 102 ) 


which satisfies the second condition in ( |100| ). The generalized-a method can be described in terms of the two 
remaining free parameters a,„ and a/. Using (|98]l and (|102|), we rewrite ^“2 in (101 1 as 


A 00 _ 

1,2 — 


ocf — a„j — 1 

a/ - a,„ -h 1 


(103) 


Let p<x> denote the user-specified value of the spectral radius in the high-frequency limit. Since we require that 
^3 < ■^ 1,2 for all D, it is p = |A“ 2 l- It has been shown that for a given level of high-frequency dissipation, i.e., for 
fixed Poo, low-frequency dissipation is minimized if ^“2 = A“, i.e., if a/ = (a,„ -(- 1) /3 (dotted line in Fig. p 
It is more convenient to describe this optimal case by defining a^ and a/ in terms of poo: 


CCifi — 


2 poc - 1 

Poo 1 


«/ = 


Poo -|- 1 


It is worth to mention that for the exact solution of (931, we get: 


Pex — 1 ) — 


1 


ft) At 


(104) 


(105) 


B Bathe-method 


Applying the Bathe-method (Fig. 331 to the single degree of freedom system introduced in ( |93] ), we can write 
the amplification matrix f7]: 



► 


Figure 33: Bathe-method in time; red: trapezoidal rule, green: Euler backward rule. 


where 


^ /-144^2-pl9D‘ 
Asathe = JZ ( —96fl^ -|- 

^ \ 144 _ 19^2 


144 ^ 2 -^ 5^4 - 28^2 \ 

144 - 47^2 48 - 4^2 

144 - 5Q? 28 / 


D= (16 + Dp (9 + Dp . 


(106) 


( 107 ) 


Analyzing the eigenvalues of Asathe shows that the method is unconditionally stable and poo = 0. 























C Floating frame of reference description for a slider-crank mechanism 


The following section presents the derivation of the system matrices for the given slider-crank mechanism in 
Fig.|34l 
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Figure 34: The slider-crank mechanism with flexible rod and unilateral constraints, element j of body i. 


C.l Finite element formulation 


For the present problem, we consider the two-dimensional beam element shown in Fig. 34 We describe the 
displacement field within the element by the following polynomials in two directions X' and Y': 

w^i = ao + a\Xi , (108) 

(109) 

I for each element, we 


Wyi =02 + ajxi + <34X1 -|- flsXi 


Using the above relation for the displacement and considering nodal coordinates (Fig. 
obtain space-independent shape functions of the beam element: 


S‘J = 


0 0^00 
0 l-3^^ + 2^^ l{^-2^^ + ^^) 0 3i§2_2^3 


( 110 ) 


where § =x/l. 


C.l Floating frame of reference 

In the floating frame of reference formulation presented in this section, the configuration of each deformable 
body in the multibody system is identified by using two sets of coordinates: reference and elastic coordi¬ 
nates B4l . Reference coordinates define the location and orientation of a selected body reference {X‘, Y‘ and 6' 














in Fig. 351. Elastic coordinates describe the body deformation with respect to the body reference {ql to in 
Fig. 34 1 . The motion of the body is defined as the motion of its reference plus the motion of the material points 



Figure 35: Deformable body coordinate. 


on the body with respect to its reference (Fig. 351. We write: 


,0 — 




( 111 ) 


where uj = {ufi M/ 2 )^ is the deformation vector of element j of deformable body i, S'^ is the shape matrix 
of element j, q'j is the vector of elastic coordinates that contains the time dependent nodal values q( to q^. For 
an arbitrary body i of the system, e.g. the flexible rod in this example, we select a body reference {X',7'}, the 
location and orientation of which with respect to the global coordinate system are defined by a set of coordinates 
called reference coordinates and denoted as q‘j.. For the planar motion of deformable bodies, which is a special 
case of three-dimensional motion, the vector q\ can be written in a partitioned form as 


q[={R^ ey 


( 112 ) 


where R' is a set of Cartesian coordinates that define the location of the origin of the body reference (Fig. 351 
and 6' is a set of rotational coordinates that describe the orientation of the selected body reference (in the present 
planar case, it is a scalar value). There is no rigid body motion between the body and its coordinate system. 
The floating frame of reference formulation does not lead to a separation between the rigid body motion and 
the elastic deformation. 

In the case of a rigid body, the global position of an arbitrary point P on the rigid body can be written in planar 
analysis as: 


r’p = R' +A‘u‘p , 

where Up is the local position of point P and A' is the transformation matrix defined as 

''cos 6' —sind'^ 


A‘ = 


sin0' cos 6' 


(113) 


(114) 


For deformable bodies, the distance between two arbitrary points on the deformable body does not, in general, 
remain constant because of the relative motion between the particles forming the body. In this case, the vector 













Up can be written as: 


= Uq + uj = Uq + 


iji 


(115) 


According to Fig. 35 and what we have discussed so far, we describe the new position of an arbitrary point P' 
on the flexible body based on reference and elastic coordinates as: 


(116) 


r'l = R‘ +A‘u/, = R‘ + A' (+ S'^qj 


We summarize all the unknowns which are necessary to calculate the new position of the arbitrary point P in 
the vector q'^: 


q 


R" 

•j = [ e‘ 


(117) 


Differentiating ( |116| ) with respect to time yields 

_L A — p' a‘„‘J a‘ 


r’J = R‘ +A u J +A‘uf = R PA u^ +A‘S'Jql , (118) 

where in case of planar motion, we have A=Aq6. Then, the velocity vector can be written as 


r‘J, = (7 A‘gu^j A’S^j) 



(119) 


v^/ 


where I is identity matrix, and A‘q is the partial derivative of the transformation matrix with respect to the 
rotational coordinate d‘: 


A'e = 


— sin 6' —cos 6* 
cos 6' —sin 6* 


( 120 ) 


Equation (119) can also be written as 


= L‘Jq‘J 


( 121 ) 


where l0 = (7 AqmO A'S'^). 


C.3 Constructing the mass matrix 

In this section, we develop the kinetic energy of deformable bodies and point out the differences between the 
inertia properties of deformable bodies that undergo finite rotations and the inertia properties of both rigid and 
structural systems. Then, we explain how to assemble the total mass matrix considering two rigid bodies (crank 
and slider) and joint constraints. 

For constructing the mass matrix the following definition of the kinetic energy is used for element j of de¬ 
formable body /: 


2 Jv'J 


( 122 ) 


where p‘-> and V‘j are, respectively, the mass density and volume of the element j, P-> is the global velocity 
vector of an arbitrary point of the element. Using the expression of the velocity vector ( 119| ), we write the 
kinetic energy as 


1 • r 


JV‘i 


p‘JL‘J L'JdV’J q’J 


(123) 













where is recognized as the symmetric mass matrix of body i. It is defined as 


/ 


M'J = / = 

Jv'i 


'VJ 


{A‘gu‘jf\{l A'qu'^ A'S'^)dV^^ 


= / . P 

JV'i 


^symm. 


A'gu'J A'S'^ 

1 dV‘J 


( 124 ) 


where the orthogonality of the transformation matrix A‘ A‘ = I is used in order to simplify the sub-matrix in the 

■T ■ ~ 

lower right-hand comer. Further, Ag A‘ = I with 


7 = 


0 1 

-1 0 


(125) 


The mass matrix (124i can also be written as 


M'J = 


rriRR 

mR9 

rriRf 


me 

mf 

symm. 




(126) 


where 


— 


-if = 
-if = 


I p‘jldV‘j , 

^ie = 

( p'-AU-dV'-, 

JV'i 

Jv'i 

[ p‘jA‘S‘jdV‘j , 

^ie = 

[ p’ju‘j\^jdV‘j 

Jv‘i 

Jvij 

[ p‘ju‘-^IS‘}dV‘j, 
Jv'i 


[ p’jS^j^S’jdV^j 
Jv'j 


(127) 


Note that the two sub-matrices and which are associated, respectively, with the translational reference 


and 

and elastic coordinates, are constant. Other matrices, however, depend on the system generalized coordinates. 
The mass matrix in the case of a rigid body motion can be written as 


m; 


rigid 


= ( 


mRR niRQ 

) symm. niQQ 


(128) 


In the case of structural systems, the reference coordinates remain constant with respect to time and the mass 
matrix of the body in this case is the constant matrix . When a deformable body undergoes rigid body 

motion, the mass matrix is defined by ( |127[ ) and the sub-matrices and nigj- represent the coupling between 
the reference motion and the elastic deformation. 

In the following, we detail each sub-matrix to find a simplified version for more efficient calculation. The 
matrix can be defined as 


mi.= lp‘jidy‘- = 


IV‘i 


m‘- 0 ( 

0 m^j I 0 ’ 


(129) 


where 7 is the identity matrix and m'- is the mass of the element j of the deformable body i. We write the 
sub-matrix as 




IV'J 


IV'J 


111 o 
Mq + “/ 


dV'J =A\ 


fi+rci} 


(130) 


where the matrices !'( and S'^ are defined as 


ri= p'-u'^dV'j , S‘^ = p‘JS'jdV‘j . 


IV'i 


IV‘i 


(131) 











The vector is the moment of mass of the body about the axes of the body reference in the undeformed state. 
It may vanish if the origin of the body reference is initially attached to the body center of mass. The vector 
qj represents the change in the moment of mass due to the deformation. Using ( |127| ), we verify that 


m 


Rf 


(132) 


The expression for m^g is 


m 


- (ij 


= m 




u I u 

+Uy 

+ ( m'Jg 


•1 T r 


O I o 

Mg +“/ 


rf 


+ 1 "^00 


dv‘j = 


ff 


IV'I 


• •T • < • • 7 ’ • < • •T > • 

Uq Uq+2uq Ujr+uj uj 


dV'^ 


(133) 


in which the sub-matrix mla reduces to a scalar that can be written as the sum of three components. The first 


component, I m 


‘00 


‘00 

is the mass moment of inertia in the undeformed state: 


m 


00 


= J p'-’uq UgdV^^ = 




,0 


dV'^ = 


(134) 


Clearly, this integral has a constant value and does not depend on the body deformation. The last two scalar 


components. 


m 


00 


rf 


and (nigg 


ff 


represent the change in the mass moment of inertia of the body due to 


deformation. These two components are evaluated according to 


m 


00 


rf Jv‘i 


= 2 I u'fdy^^ = 2 


m 


ee 


= p’ji/fi/jdV‘J = q^f 


ff JVJ 


IV'I 


JVJ 


p’jufs'jdV'J 


O IJ 


qj. = 2 qj 


p'^S‘J‘ S‘JdVJ q'j 


(135) 

(136) 


If we use definition (127), we can write the following: 


O \ O O 


(137) 


where 




(138) 


Finally, we introduce 


<f = jvf 


I o 

Uq -\-uj 




(139) 


where the constant skew symmetric matrix S ^ is defined as 


I'J = / p‘Ju‘J I S’^dV^j , ~S'^= p‘js’j I s’jdv‘j. 
Jv‘j Jv‘j 


(140) 


We conclude that, to completely describe the inertia properties of the deformable body in plane motion, a set 
of inertia shape integrals is required. These integrals, which depend on the assumed displacement field, can 
be obtained using the Gaussian quadrature method. As we are dealing with polynomials for describing the 
displacement field, for finding the optimum number of Gaussian points, we have to know the highest degree of 
the polynomials which appear in the calculation of the mass matrix. We expect polynomials with degree 6 at 
most, which need 4 Gaussian points to yield an exact integration. Figure shows Gaussian points which are 
necessary to evaluate exact values on a sample element. To evaluate the integrals for § G [0,1], we have to map 
Gaussian points and weights. 

Once we have calculated the mass matrix on element level, we have to assemble the total mass matrix consid¬ 
ering mutual degrees of freedom. One should note that entities in the mass matrix of different elements share 
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Figure 36: Gaussian point distribution for exact integration of polynomial representations in mass matrix eval¬ 
uations. 


the same reference coordinates, but perhaps different elastic coordinates. Figurej^shows the regions in which 
we have the overlap (purple color) between element j (light red) and element j +l (light blue) in the assembly 
process. As the slider and the rod are attached to each other at one revolute joint, the end point of the flexible 
rod and the mass center of the slider have the same translational velocity. Therefore, we add the effect of the 
rigid slider to the mass matrix of the last element of the flexible rod, also by adding one additional degree of 


freedom in the total mass matrix for 63 (Fig. 341: 


q‘’^ 






(141) 


where index 3 is related to the slider (third body), and I is related to the last element of the mesh which is 


connected to the slider by means of a revolute joint. We assemble the first term of ( |141| ) in the total mass matrix 
regarding mutual degrees of freedom and we add one additional row and column for the new degree of freedom 
63 which contain zero everywhere except one diagonal term which is J 3 . 

To add the effect of the rigid crank mass into the total mass matrix, first we define the constraint for the 
connecting joint between crank and flexible rod: 


= /?' - 


Zi cos 6 
11 sin 6 ] 


= 0 , 


(142) 


where = (+' y') is the translational coordinate of the rod reference frame. Equation ( |142 ) shows that 
the rod reference coordinate and therefore its derivative with respect to time can be expressed in terms of 61 . 
Therefore and for rewriting the total mass matrix in terms of 61 , we need the time derivative of the constraint 





























Figure 37: Mass matrix results from assembling of element mass matrices j and 7 + 1; purple: summation of 
mutual DOFs, red: element mass matrix j, blue: element mass matrix 7 + 1. 


C': 


c' = r ‘ - 


/—Zi sin 0 i\ 

\ ZicosFi J 


di =R‘-Cl^ei=0. 


(143) 


Using ( |143| ) and keeping in mind the kinetic energy formulation ( |123| ), 
that it only depends on the coordinate di according to the constraint : 


we modify the total mass matrix such 


+/i , 

mRo = Cqj niRQ , theR = , (144) 

ifiRf = Cg^niRf , fhfR = fhlf , 


where I\ = Ji + ^m\l^ is the rigid crank mass moment of inertia with respect to the joint. We consider the 
contribution of the rigid crank in the total mass matrix according to the kinetic energy of the crank in terms of 
its only degree of freedom 61 : 


T^ = \heh (145) 

where index 1 denotes the body number of the crank. 

Figure]^ shows the final mass matrix configuration after considering all the effects in the slider-crank mecha¬ 
nism. 
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Figure 38: Total mass matrix for the slider-crank mechanism considering joint constraints. 

C.4 Quadratic velocities 

For the quadratic velocity vector, we need to calculate the derivative of the mass matrix with respect to time 
and the derivative of the kinetic energy with respect to the degrees of freedom. For the time-derivative of the 
mass matrix, we write 


= I p‘J ( V^‘ L‘J) dV'J = 


rriRR 

mRQ 

niRf 


rhee 

thgf 

symm. 


'«//, 


(146) 


where and are the modified versions of and according to the hinge constraint defined in ( |142| ): 


Dj = cL e, di A'eS'^cij - A' u‘^ + 62 A'^S^^dz 


(147) 

(148) 


Thereby, Cg g is the second derivative of the joint constraint: 


ce. = (7.:r:o:)®i=CM.ei 


Zi cos 61 
11 sin 61 


(149) 


After assembling the element contributions according to the pattern described in Fig. we add the effect of 
the rigid slider to the matrix of the last element of the flexible rod, noting that J 3 is zero: 


,-.3 

M =m^ (L L +L L 


( 150 ) 


Note that the effect of the rigid crank is already taken into account by modifying the L matrix. 





















































































To calculate the second term of the quadratic velocity, first we write the expression for the total kinetic energy 
as 


T = ^mRRQl + 0imRed2 + + d\mRfqj: + d2tnQfqf + -ntffqj: + -Jsd^ 


(151) 


Keeping in mind definition (144 1 , the derivative of the kinetic energy with respect to the nodal coordinates is 


dT 

dq 


(cl" {m 2 + m 3 ) Cl^ 01 + (/i + %.) 02 + 

-C^"a (/i +Sqj) eid2+c^J^Alseiqf 

0 

C^"A'e50i02 + (/3 +%)^/0| + 5024/ 


where S, S//, S, /i, and Ii, are the assembled versions of S‘L, S \ l'( and defined in 


C.3 


(152) 


Finally, we have 


Qv = - 


/ mRR 

mRQ 0 

mRf\ 

/ 0 i\ 


ihge 0 

the/ 

02 


0 

0 

03 

\symm. 


'«/// 

\9// 


+ 


dT 

dq 


(153) 


C.5 Stiffness matrix 

Considering a linear isotropic material, the virtual work due to the elastic forces for element j of body i can be 
written as 

5wi^ = - [ o‘j^de'jdV‘j, (154) 

Jv'j 

where a'^ and e'^ are, respectively, the stress and strain tensors. Since the rigid body motion corresponds to 
the case of constant strains and since we defined the deformation of flexible bodies with respect to the body 
reference, the strain displacement relations can be written in the following form: 

e‘j=D‘ju‘i, (155) 

where D‘^ is a differential operator. We write 

e‘j = D’js‘Jq‘i . (156) 

For a linear isotropic material, the constitutive equations can be written as 

= (157) 


where C‘^ is the symmetric matrix of elastic coefficients. We conclude 

= C^D‘js‘jq‘i , (158) 

where the stress tensor is written in terms of the elastic generalized coordinates of body i. Finally, we have 


5Wl^ = -q} 


IV'i 




5q'l = -q‘l K'l^.5q‘} , 


where 




(159) 


(160) 


Neglecting the shear deformation and using the assumptions of Euler Bernoulli beam theory, the strain energy 
for the element j of the elastic rod can be written as 




rl'J 


(U/j ^f2} 


ij fE'iA'i 

( 0 


0 

Edid 


cu 


J \ ‘J 
/I 




dr 


(161) 














where l'-> is the length of element j of the beam, E‘J is the modulus of elasticity of the beam, A‘-> is the cross- 


sectional area, is the second moment of area, 


and ^<^2 are the axial and transverse displacements of 


element j respectively and (') denotes differentiation with respect to the spatial coordinate. Taking the derivative 
of ( 161| ) with respect to < 7 ^ and comparing with ( 159| ), we write 


2 ^ 0 

\ 0 E^jpj 


Using the shape functions introduced in (llOi, we conclude: 


DiJsU = s'‘j = I g ^ 

V ^ i-P- 



0 0 10 0 

-6 + 12^ pj{-4 + 6^) 0 6 - 12 ,^ 


(162) 


(163) 


where the first matrix is the inverse of the Jacobian for transforming to local coordinates with ^ =x/l, and the 
second matrix is the derivative of 5"^ with respect to ^. We calculate the stiffness matrix for each element using 
( |160| ). It is worth to mention that as we are dealing with derivatives of the shape functions, the highest degree 
of polynomials which appears in the integral of ( |160| ) is 2. Therefore, we only need 2 Gaussian points for the 
computation of the stiffness matrix on element level. Figure shows the stiffness matrices for element j and 
7 -|- 1 and the assembling procedure. Blank areas, which correspond to the reference degrees of freedom, i.e., 
61 , 02 and 63, are filled wifh zeros. 



Figure 39: Stiffness mafrix assembling of fwo elemenf j and 7 + 1 stiffness matrices; purple: summation of 
mutual DOFs, blank: zero entries. 


C.6 External forces 

The virtual work of all external forces F' acting on element 7 of body i in the multibody system can be written 
as 

[ E‘j^5E^dV^j, 

Jv'i 


( 164 ) 








































where 


/5R'' 


5r' = {l A'gu^j A'S'i)\ 5e^\ =Vi5q'^ . 

\8ql 


( 165 ) 


Combining the above two equations, we get the following expression for the virtual work of external forces: 


where 


'5R' 

sw :^={ q ^ ef ef) 1 

Ml 


ef = / Mdv^j , 

JV‘j 

Q‘f = [ A'qu'^AV'^ 


(166) 


'-e — I ^ 

JV'J 

Qf = / 


(167) 


IV'} 


Assembling of the external force vector can be done with the same procedure as for the mass matrix. To add 
the virtual work of the rigid slider, we use the revolute joint constraint between slider and rod 


Qf =F^V\ 


(168) 


where index 3 is related to the slider (third body), and index I is related to the last element of the mesh which 


is connected to the slider by means of a revolute joint. We assemble ( |168[ ) in the total external force vector 
regarding mutual degrees of freedom and adding one additional row and column for the degree of freedom 63 
which is zero as there is no moment in this direction. To add the effect of the rigid crank mass into the total 


external force vector, we use the constraint for the connecting joint between crank and flexible rod (142). We 


modify the sub-vector introduced in (167 1 , such that it only depends on coordinate 61 according to the constraint 
C': 


Qf =Qfcl^-mgMne,, (169) 

where the second term comes from the work of the gravitational force with the rigid crank. It is worth to 
mention that according to the polynomial order in the external force integral, 2 Gaussian points are sufficient 
for each element. 


C.7 Unilateral constraints 


Figure shows the geometric characteristics of the translational clearance joint. It is 2a the length and 2b 
the height of the slider. The height of the notch is given by d. The existence of a clearance in a translational 
joint introduces two extra degrees of freedom. Hence, the slider can move freely inside the guiding limits until 


it reaches the surfaces. Considering the geometry of the slider according to Fig. 41 and keeping in mind the 
revolute joint between slider and end point of the rod, we write the position of the slider center according to 


(116 1 and (142i as 


’ eg — ( / ) —F +A Ucg — 

.teg 


h cos 01 
Zi sin 01 


+ A‘(02) 




(170) 









Figure 40: Different scenarios for the slider and notch interaction. Red points indicate active contact points. 



Figure 41: Definition of the gap functions for the slider-crank mechanism with unilateral constraints. 

where index I is related to the last element of the mesh of the flexible rod wifh initial lengfh I 2 . Therefore, fhe 
nonlinear normal and fangenfial gap funclions splif up for each comer: 


SNi {q) = ^- rig + a sin 63 - f 7 COS 63 , 

(171) 

8 N 2 {q) = ^-r^g-asine3-f7cos03 , 

(172) 

SNiiq) = ‘^+rlg-a&me-i-bcosdi, , 

(173) 

gN^q) = ‘^+rlg+a&me-i-bcosd 2 , 

(174) 

gTi {q) = rig -acosd 3 -b sin 63 , 

(175) 

gT 2 (<?) = rlg + a cos 63 -b sin 63 , 

(176) 

gTi (<?) = rig - a cos 03 + sin 63 , 

(177) 

gT^ {q) = rig + acosd 3 + b sin 63 . 

(178) 




























The matrices of generalized force directions are the derivatives of gj^ and gj with respect to q: 


K{q) 


Wl{q) 


/ — l\ cos 01 

— [AgMcg]^ 

a COS 03 -|-Z?sin 03 


—Zi cos 01 

— [AgMcgjy 

—a cos d^ + b sin 03 

-A 4, 

Zi cos 01 

\AQUCg\y 

—a cos d^ + b sin 03 


y Zicos0i 

[AgMcg] 

a cos 03 -|-Zjsin 03 

A4, / 

/—Zi sin 01 

[AgMcg]^ 

<2 sin 03 — Z? cos 03 



—Zi sin 01 


-a sin 03 — Z? cos 03 



—Zi sin 01 

\A^qUc^^ 

a sin B^ + b cos 03 



\—Zi sin 01 

[AgMcg]^ 

-a sin 03 -I-Z 7 COS 03 




( 179 ) 


(180) 


where [= 1 =]^ and [=t=]^ mean the component of the vector respectively in x— and y—direction. The matrix S is 
defined as follows: 


S = 


0---0 0 0 0 1 0 0 \ 

0---0 0 0 0 0 1 0 „ 

/ IxN 


(181) 


where N = 3 x {neu +1) is the number of elastic degrees of freedom. The last six columns are filled with 
= 1) related to the degrees of freedom of the last element of the mesh. 


C.8 Boundary conditions 


When solving a finite element problem with the floating frame of reference formulation, the system becomes 
singular. This is because rigid body motion is added to the equations of motion at the same time as the defor¬ 
mation field also contains rigid body motion. A modal analysis would show that the first three eigenvalues are 
equal to zero which correspond to the degrees of freedom of the rigid body in planar motion. In this model, 
boundary conditions for the body reference system are introduced to avoid the singularity in the system of equa¬ 
tions. There are different ways to define the boundary conditions. Shabana Il43]| shows that two sets of modes 
associated with two sets of boundary conditions can be used to obtain the same solution if the coordinate system 
is properly selected. Therefore, the physical deformation is unique in the inertial frame. For a beam, it is most 
common to use a clamped-free or a simply supported reference system (Fig. 42 1 , where three conditions are 
given in both cases. Clampled-free, i.e., tangential, means that the reference system is tangential to the beam 



Figure 42: Tangential and pinned reference system. 


deflection at the root of the beam, i.e., both displacements and rotations are equal to zero: 

9/1 ~ 9/2 = 9/3 ~ 0 • ( 182 ) 

For a simply supported, i.e., pinned, reference system, the root of the beam is locked and the end of the beam 
is moving but only along the local x direction: 

9/1 =9/2 = 9/5=0, (183) 

where index I stands for the last beam element. In the present simulation, the tangential reference system is 
chosen. We compare different boundary conditions using modal coordinates in Sect.|4.5|as introduced in[D| 














D Modal analysis 


Modal analysis is a process where the nodal displacement vector is approximated by a linear combination of 
dominant eigenvectors (also called mode shapes) as it is shown in Fig. 43 Elimination of high-frequency mode 
shapes decreases the number of numerical operations per time step because the size of the matrices in the 
equations of motion is much less than in the non-reduced case. 
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Figure 43: Superposition of the modes in modal analysis. 


For the flexible beam with mass matrix niff and stiffness matrix Kff, we obtain the angular frequency (Ok and 
the relative mode shape (/»using the free vibration equations of motion. The free vibration equations of motion 
can be derived from ( fT^ in case of no external forces and damping: 

{Kff-(oj niff) = ^ ■ ( 184 ) 

The deformations are expressed as 

nm 

= (185) 

k=i 

where qf - is the deformation of the degree of freedom number j. Considering that the transformation matrix <I> 
consists of the eigenvectors (j) the values ^jk and are components of the eigenvectors and modal coordi¬ 
nates (the new unknowns). The number of reduced coordinates is chosen depending on the accuracy. We 
rewrite the above equation: 


qf = <!><?„ 


The deformation field (|111|) is described with 


7^ - 






(186) 


(187) 


Using characteristics of normalized orthogonal eigenvectors ((/»,(-), we decouple the sub-matrices niff and Kff 
in ( 160|) and ( 138[) respectively: 


^^niff^ = ln„ 


x-nm ; 


^‘Kff^ = (ot 8 


ki 


(188) 


In order to complete the decoupling in the mass matrix, we deal with the sub-matrices niRf and rngf on 
element level using the modified shape functions introduced in ( 187| ). Moving the body coordinate system to 
the mass center, we show that the integrals in (|131|) and (|140|) vanish if 


f (l)k{r)dV =0 , [ r^k {r) dU = 0 . 

Jv Jv 


(189) 


It can be shown that the free-free modes satisfy the condition above and also the mean-axis conditions which is 
obtained by minimizing the kinetic energy of the elastic motion with respect to an observer sitting on the flexible 
body |[3. This ideal coordinate system and the free-free eigenfunctions are shown in Fig. 44 Removing the 
first three modal coordinates (corresponding to the rigid body motion) results in decoupled versions of the total 
mass matrix M. It is clear from Fig. that the deformation at the two ends of the beam do not vanish as 
defined in this coordinate system. Shabana 1431 modifies the free-free shape functions to satisfy the boundary 
condition which results in a simply supported, i.e., pinned, mode. However, it cannot be used to decouple the 
mass matrix anymore. Using articulated-free modes as shown in Fig. [45]l[T^ . i.e., fixing the position of the left 
end and leaving the right end free, we can fix the boundary condition at the joint in addition to satisfying the 
second condition in (|189|), which results in vanishing terms regarding nigf. 































Figure 44: Free-free boundary condition and new reference frame at the mass center for decoupling the refer¬ 
ence and the elastic coordinates. 



Figure 45: Articulated-free boundary condition for partially decoupling the reference and the elastic coordi¬ 
nates. 
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